Spontaneous traveling waves in oscillatory systems with cross diffusion 
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We identify a new type of pattern formation in spatially distributed active systems. We simulate 
one- dimensional two-component systems with predator-prey local interaction and pursuit-evasion 
taxis between the components. In a sufficiently large domain, spatially uniform oscillations in such 
systems are unstable with respect to small perturbations. This instability, through a transient 
regime appearing as spontanous focal sources, leads to establishment of periodic traveling waves. 
The traveling waves regime is established even if boundary conditions do not favor such solutions. 
The stable wavelength are within a range bounded both from above and from below, and this range 
does not coincide with instability bands of the spatially uniform oscillations. 

PACS numbers: 87.10.+e, 02.90.+p 



I. INTRODUCTION 

Dissipative structures, i.e. patterns in spatially ex- 
tended systems away from equilibrium have been inten- 
sively studied for many decades. A very comprehensive 
review can be found in Cross and Hohenberg [1]; results 
obtained since then would probably require an even more 
extensive review. A very popular class of mathemati- 
cal models is the reaction-diffusion systems with diag- 
onal diffusion matrices. There have been numerous in- 
dications that non-diagonal elements in diffusion matri- 
ces, i.e. cross-diffusion, can lead to new nontrivial effects 
not observed in classical reaction-diffusion systems, e.g. 
quasi- solit oris in systems with excitable reaction part [2- 
6] . However oscillatory systems are more prevalent than 
excitable and nontrivial effects of cross-diffusion in oscil- 
latory systems have not been studied yet. Here we con- 
sider an example where the reaction part of the system 
is dissipative while the diffusion part is not. We describe 
spontaneously generated periodic waves, and identify the 
features of these waves that indicate that we are dealing 
here with a phenomenon not seen before. 

A general formulation of a reaction-diffusion system 
with nonlinear diffusion is 



V(D(u)Vu), u,f G 



pN 



D G 



pNxN 



Both the reaction term f(u) and the diffusion term 
V(D(u)Vu) in the right-hand side represent dissipa- 
tive processes, which for diffusion implies that matrix 
D G R 7VxAr is positive (semi-)definite, typically diagonal 
with non-negative elements. A huge amount of results 
have been obtained about pattern formation described 
by such models. However, many physical situations lead 
to non-diagonal elements in D, i.e. cross-diffusion (see 
e.g. discussions in [7, 8]). Some such situations may be 
adequately described by D whose eigenvalues have zero 
real part, e.g. when the self-diffusion of components is 
negligible. In such cases reaction part is dissipative and 



the "diffusion" part is not. Physical consequences of such 
ambivalence are little understood yet. 



Cross diffusion has been seen to produce interesting 
phenomena, such as fronts, pulses and stationary periodic 
structures (see e.g. [9, 10] among many other works), 
however phenomenologically similar regimes are known 
in reaction-diffusion systems, too. 



In a recent series of works we have described unusual 
phenomena, such as quasi-solitons and their variations, 
in excitable systems in which linear or nonlinear cross- 
diffusion was added to or replaced self-diffusion (see e.g. 
[2-6]). The ability of a medium to conduct solitary waves 
is stipulated by its excitable kinetics described by the re- 
action term f (u), whereas specifics of their interaction are 
also due to the cross-diffusion terms. However, excitabil- 
ity is a relatively exotic, albeit very important, type of 
behaviour compared to oscillations. For instance, in pop- 
ulation dynamics, plausible excitable predator-prey mod- 
els have been proposed [11] but we are not aware of reli- 
able observations of natural systems described by such 
models. On the other hand, oscillatory behaviour in 
predator-prey systems is textbook material [12, 13] and 
there are plentiful observational data on traveling waves 
in cyclic populations [14]. 



Solitary waves in oscillatory systems are not feasible, 
and it is not clear what new features cross-diffusion may 
impose. 



The purpose of this communication is to describe 
new phenomena we have observed in oscillatory systems 
with "pursuit-evasion" nonlinear cross-diffusion interac- 
tion between the components. 
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FIG. 1: Different regimes resulting from random perturbaton of uniform oscillations. Shown are density plots: space x is 
horizontal, time £ is vertical increasing upwards, u = 1 corresponds to black and u = corresponds to white, (a) TB model, 
L — 15, taxis (h- = 1, /i+ = 0, .D n = D v = 0), periodic boundary conditions, £ G [0,2500]. (b) Same, no-flux boundary 
conditions, £ G [0,1200] U [43800,45000]. (c) Same as (a) except D u = D v = 0.05, h+ = h- = 0. (d) Same as (a) except 
w = 0.07, L = 10, £ G [0, 1200] U [36300, 37500]. (e) Same as (a) except h+ = 0.1, L = 50, £ G [0, 1200] U [5000, 6200]. (f) Same 
as (a) except h+ = 0.1, D u = D v = 0.02, L = 50, £ G [0,400] U [1900,3900]. (g) RM model, /i_ = I, h+ = D u = D v = 0, 
L = 25, £ G [0,2500]. 



II. THE MODELS 

We consider two predator-prey models with cross- 
diffusion terms of "pursuit-evasion" mutual taxis, 

du p/ x _ <9 2 w , d f dv\ 



dv ) + D ^ V h ^ ( ^ U 

dt ' v <9x 2 dx V <9x 



(2) 



for (x,i) G [0, L] x [0, £ ma x] for two reaction kinetics, the 
Truscott-Brindley (TB) model [11] 



f(u,v) = -u) -vu 2 /(u 2 + ^ 2 ) ? 

g(u,v) = jvu 2 /(u 2 + v 2 ) — wv, 



(3) 



where /? 0.43, 0.053, 7 0.1 and w 0.055 unless 
stated otherwise, and the Rosenzweig-MacArthur (MA) 
model [13-15] 



f(u,v) = /3u(l — u) — vu/(u + z/), 
g(u,v) = jvu / (u + v) — wv , 



(4) 



where /? = 1, z/ = 0.3, 7 = 0.15 and u> = 0.03 unless 
stated otherwise. Here u represents prey, v predators, 



the term with h+ describes pursuit of prey by predators 
and the term with h- describes evasion of predators by 
prey. The simulation were done on an interval x G [0, L] 
with periodic or Neumann boundary conditions for both 
components, using forward Euler stepping in time, center 
differences for the diffusion terms and upwind difference 
for the taxis terms, see Tsyganov et al. [3] for details and 
justification. Except where stated otherwise, we used 
discretization steps Ax = 0.1 and At = 4 • 10 -4 . 



III. NUMERICAL OBSERVATIONS 

Figs. 1 and 2 illustrate the phenomenon of the spon- 
taneous onset of periodic waves. Starting from arbitrary 
spatially-uniform intial conditions at £ = 0, after a tran- 
sient allowed to establish uniform oscillations, pertur- 
bations were introduced and subsequent evolution ob- 
served. The perturbation was introduced at half of the 
grid points chosen randomly, where at £ = 300 the values 
of u were replaced by randomly chosen numbers in the 
interval between 0.15 and 0.45. Fig. 1 shows space-time 
density plots and fig. 2 illustrate selected profiles of the 
emerging wavetrains. 

In the TB model with periodic boundary conditions 
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FIG. 2: (color online) Examples of profiles of spontaneously 
established periodic waves. Shown are dependencies of u and 
v on x at fixed £, and direction of propagation by arrows. 
Parameters are the same as in fig. 1 except interval lentgh 
L, specifically, (a) as in fig. 1(a), (b) as in fig. 1(d), (c) as in 
fig. 1(e), (d) as in fig. 1(f) and (e) as in fig. 1(g). 

(fig. 1(a)), after a "random" transient lasting two or three 
bulk oscillation periods, patterns start to emerge: waves 
start "from nowhere" and annihilate upon collision with 
other such waves. After a few periods of such collisions, 
the waves propagating leftwards win over and a periodic 
wavetrain establishes which then persists. Different seeds 
in the random number generator produce solutions dif- 
fering in details but always leading to periodic trains, 
leftward and rightward propagating with equal probabil- 
ity (compare density plots and wave profiles in fig. 1 and 
fig. 2, which corresponded to different simulations with 
the same parameter sets). 

Impenetrable boundaries do not allow periodic wave- 
train solutions; however the tendency to establish peri- 
odic wavetrains is observed even then. In fig. 1(b) right- 
ward propagating waves win over. Their impact with the 
right boundary x = L is with partial reflection, when the 
reflected wave is weak and soon decays; note that this be- 
haviour is typical for collision of solitary excitation waves 
in such systems [3]. The left boundary has a quenching 
effect, but at a distance from it waves emerge sponta- 



neously. This distance varies irregularly, indicating that 
spontaneous generation of waves is associated with an in- 
stability, thus sensitive dependence on initial conditions 
and probably chaotic dynamics. This irregular pattern 
persists for a long time. 

This behaviour is in a contrast with a system with 
the same kinetics but pure diffusional spatial terms: in 
fig. 1(c), similar initial random perturbations lead very 
quickly to re-establishment of spatially uniform oscilla- 
tions. 

The parameters used in fig. 1(a) are close to the bound- 
ary of the oscillatory regime in the TB model (achieved 
e.g. at w ~ 0.053 with other parameters fixed). When 
parameters are further into the oscillatory region, sponta- 
neous generation of periodic wavetrains is still observed, 
although the transient period of spontaneous wavelet 
generations and collisions lasts longer, see fig. 1(d). 

We have also found that prevalence of the "evasion" 
taxis (h- coefficient) helps generation of periodic trains, 
but h+ = is not necessary, and such generation can 
be observed with the "pursuit" taxis present as well, see 
fig. 1(e). 

Spontaneous generation of periodic trains is observed 
in the RM model as well, see fig. 1(f). 

The spontaneously emerging periodic wavetrains typ- 
ically had wavelengths in a limited range. To check 
whether this is dictated by initial conditions or is due to 
limitations of the system, we performed simulations in a 
circle, i.e. an interval with periodic boundary conditions, 
of a slowly changing length L. We started from an estab- 
lished propagating wave in a circle. Then we changed the 
length L of the circle in small steps, allowing sufficient 
time between the steps for the waves to adjust. During 
the simulation we monitored the number of waves n, de- 
termined via the number of points where u crossed the 
level u = u* = 0.2, and the periods T defined as inter- 
vals between u crossing the level u = u* in the positive 
direction. Results of one such simulation are shown in 
fig- 3. 

The number of waves n in the interval did not remain 
constant (fig. 3(a)), but spontaneously adjusted so as to 
keep the average wavelength within certain limits: be- 
tween approximately 2.5 and 8 in the simulation shown. 
This number was not a unique function of the interval 
length: changing L upwards and downwards produced 
different dependencies n(L), i.e. we have hysteresis. Sim- 
ulations at slower rate of change of L slightly changed the 
n(L) dependencies but the hysteresis stayed. Near the 
transition points where n changed the value, the propa- 
gation of the waves was non-stationary, and was always 
for L just below the transitional value, whether it was de- 
creasing (fig. 3(b)) or increasing (fig. 3(c)). Increasing L 
had a noticeably more destabilizing effect than decreas- 
ing. 

The nature of the non-stationary solutions is illus- 
trated by the density plots shown in fig. 3(d). Starting 
from an n = 1 solution, an increase of L above the value 
of L w 5 leads to an instability of the steady propagating 
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FIG. 3: (color online) Variability of the wavetrains at chang- 
ing L. (a) Number of waves in the interval [0, L] as L grad- 
ually increases (by 0.2 every 2000 time units, solid red line) 
and decreases (by 0.2 every 1000 time units, dashed blue line). 
Oblique dashed lines: n = L/2.5 and n — L/8 to guide the 
eye. (b) Wave periods measured at a point, as function of L as 
it decreases, (c) Same, as L increases, (d) Density plots of two 
episodes of simulation of 1200 time units duration each, with 
L increasing by 0.2 every 1000 time units. Lower episode: soft 
transition from steady 1-wave solution to modulated 1-wave 
solution (L : 4.8 — > 5.6). Upper episode: subsequent sudden 
transition from modulated 1-wave solution to a steady 2- wave 
solution (L : 7.2 -> 7.6). 



wave solution. This is a soft, Eckhaus-type instability 
and leads to a mild modulation of the wave, producing 
a seemingly two-periodic motion. The amplitude of the 
modulations grows as L increases, until at L = 7.6 a qual- 
itative transformation occurs. A gap between the wave 
and its own copy around the circle grows so big that at a 
certain moment it is sufficient to allow spontaneous gen- 
eration of another wave, leading to an n = 2 solution. 
This solution is steady, i.e. propagates without modula- 
tions, until L grows so big it in turn becomes unstable 
etc. 



IV. PRELIMINARY THEORETICAL 
CONSIDERATIONS 

Substantial theoretical analysis of the phenomenon of 
the spontaneous traveling periodic waves is beyond the 
scope of this communication. Here we consider one naive 
approach and then some known pattern formation mech- 



anisms, which a priori might look relevant to this phe- 
nomenon, only to eliminate them, as not providing a sat- 
isfactory explanation. We will refer to the historical re- 
view by Cross and Hohenberg [1, p. 870] (CH for brevity), 
and to a recent symmetry based classification of instabil- 
ities and bifurcations of periodic dissipative waves and 
structures given by Rademacher and Scheel [16, p. 2680] 
(RS for brevity). 

a. It is not captured by lambda- omega approach The 
simple class of two-component reaction-diffusion sys- 
tems introduced by Kopell and Howard [17] and called 
"lambda-omega systems" , and closely related to the com- 
plex Ginzburg-Landau equation, allows exact solutions 
in the form of periodic waves. It has offered qualitative 
insight in many nonlinear wave phenomena, including pe- 
riodic waves in cyclic populations [14]. However, it does 
not seem to be helpful in our present case. The modifica- 
tion of the lambda-omega system, corresponding to the 
choice of signs of taxis terms in (2) is 



dz 
di 



(A(\z\) + iQ(\z\))z-i\7 2 z 



(5) 



where z is a complex dynamic variable representing u + 
iv, and the purely imaginary diffusivity here corresponds 
to the absence of self-diffusion, D u = D v = 0. Then 
the periodic traveling wave ansatz z = a exp[z(o;t — kx)], 
a,uu, k G M, gives the finite system 

A (a) =0, uu = ft (a) - /c 2 , 

i. e. all waves have the same amplitude which is a root of 
A(), and exist for all wavelengths k rather than in a finite 
interval. Stability analysis and consideration of nonzero 
self-diffusion do not help either. 

b. It does not emerge via Turing mechanism. The 
instability of spatially-uniform solutions in favour of non- 
oscillatory, spatially-periodic solutions with periods in a 
finite range is, of course, a defining feature of the Turing 
patterns, called just so by RS and classified as type I s 
in CH nomenclature. Cross-diffusion can provide an al- 
ternative to the original Turing's short range inhibition - 
long range activation condition. Indeed Turing- type in- 
stabilities and spontaneously occurring, self-supporting 
time-stationary spatially periodic patterns have been ob- 
served in locally multistable systems with cross-diffusion 
[9]. Our present observations are different in that here 
we are dealing with time-oscillating phenomena not just 
space-oscillating. 

c. It does not emerge via Turing-Hopf mechanism. 
A Hopf bifurcation of the spatially uniform equilibrium 
at a nonzero wavelength, is called "Hopf", "oscillatory 
Turing" and "Turing Hopf" instability by RS, classifed 
as type I Q in CH nomenclature, and also known as 
short-wave instability or finite-wavelength instability. It 
can lead to stable periodic propagating waves, in lasers, 
fluid convection and reaction-diffusion models [18-21]. 
In reaction-diffusion context, such waves have been ob- 
served experimentally and in simulations in populations 
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and BZ reaction [22, 23]. However, the standard way 
such instability occurs in systems (1) implies existence 
of an equilibrium that is stable with respect to spatially 
uniform perturbations, which we do not have here, and 
it only can occur if N > 3 whereas we have only two 
components, u and v. 

Specifically, for u(x, t) = u r + v e xt+lkx where u r = 
(u r ,v r ) is the spatially uniform equilibrium and |v| <C 1, 
we have the characteristic equation 

det (F r - T) r k 2 - Al) = 0, 



where F r = F(u r ) = (df/du) u= 



is the 



fn fl2 
/21 Hi 

Jacobian matrix of the reaction terms and D r = D(u r ) is 
the diffusion matrix, both evaluated at the equilibrium. 
Considering for simplicity the cases of fig. l(a,b,d,f) 
h-U r 



where D r 







we have 



A = \ (/ll + f 2 2 ± Vifll ~ f22? + 4/12/21 - /2lh-Urfc 2 ) 

which for k 2 > max (((/ n - / 22 ) 2 + 4f 12 f 2 i)/(h-U r ) , 0) 
gives oscillatory behaviour of perturbations, but then 
Re(A) = (/11 + /22V2 = const whereas it has to have 
a maximum at a positive k 2 for this mechanism to be 
relevant. 

d. It does not emerge via Turing-Hopf instability of 
spatially uniform oscillations The next possible candi- 
date is the instability of spatially uniform oscillations 
with respect to perturbations which nonzero freqency 
and nonzero wave number. This case is not considered 
in the CH nomenclature, and is called "Hopf" instabil- 
ity of spatially homogeneous oscillations, with the same 
variants as in the previous case, by RS. This instability 
looks plausible as spatially homogeneous (spatially uni- 
form) oscillations in our systems are indeed possible and 
even stable in small spatial domains, so we have investi- 
gated this possibility with particular care. As limit cycles 
in the point systems of (3) and (4) can not be described 
analytically, the investigation of stability has to be done 
numerically. We have considered solutions of the form 
u(x,t) = u Q (t) +Re (y(t)e ikx ) with |v| <C 1, which gives 
a coupled system of ordinary differential equations 



du 

~df = f(Uo) ' 

^ = (F( Uo ) - D(u )fc 2 ) v, 



(6a) 
(6b) 



with parameter k. We solved system (6) forward in time 
with initial conditions for bulk oscillations u in the basin 
of attraction of the limit cycle, and arbitrary nonzero 
initial conditions for the perturbation v. Then we es- 
timated the Lyapunov exponent for the v-subsystem, 
A(fc) = lim^oo t _1 ln(| |v(t)| |). The estimation was done 
by finding maxima of the first component of v(t) and lin- 
early fitting their logarithms against t, for an interval of 
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FIG. 4: (color online) Emergence of spontaneous periodic 
waves through instabilities, (a) Spectrum X(k) of harmonic 
perturbations of the spatially uniform oscillations (black lines 
with points) , and the histograms of the wavenumbers of spon- 
taneous wavetrains in the simulations shown in fig. 3 for in- 
creasing L (red solid line) and decreasing L (blue dashed line), 
(b) Same, for the RM model, parameters as in fig. 1(f), his- 
tograms obtained by increasing L by 0.2 every 2000 time units 
from 5 to 50 and decreasing it back to 2 by 0.2 every 1000 
time units, (c) Emergence of standing periodic waves via an 
instability of the spatially uniform oscillations. Parameters as 
in fig. 1(a) and fig. 3, L = 12.26, Ax = L/63, At = 5 • 10" 5 , 
t G [0,1000]. (d) Emergence of spontaneous periodic wave- 
trains via an instability of periodic standing waves. Continu- 
ation of (c), t e [44460, 45460]. 



large enough values of t. For selected values of k we used 
two linear independent sets of initial conditions for v, to 
eliminate the theoretical possibility of accidentally chos- 
ing initial conditions that did not lead to the maximal 
exponent. 

The resulting graph A(fc) for the TB model at the same 
parameters as in fig. 1(a) and 3 is shown on fig. 4(a). 
For comparison, we also show histograms of the empiri- 
cal wavenumbers observed in simulations shown in fig. 3, 
calculated as k = 27rn/L, separately for the growing and 
decreasing L. Fig. 4(b) shows similar graphs made for the 
RM model at the same parameters as in fig. 1(f). It is 
clear that, although there are finite bands of wavenum- 
bers producing growing perturbations, the actually se- 
lected wavenumbers are not the same as those of the 
fastest growing perturbations, and for the TB model they 
even partly fall in the interval of decaying perturbations. 
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Moreover, the growing perturbations of the spatially 
uniform oscillations in fact do not represent propagating 
periodic waves, but standing waves. This is illustrated 
in fig. 4(c) where we show a density plot of a simulation 
of the full model, similar to fig. 1(a) but with different 
initial conditions. Here we chose initial conditions as spa- 
tially uniform oscillations plus a very small perturbation 
sinusoidal in space. Note that for the limit of infinitely 
small perturbation amplitudes this exactly corresponds 
to system (6). 

We conclude that although the cross-diffusion driven 
instability does indeed take place in the considered ex- 
amples, the waves that emerge are in fact quite different 
from the spontaneous periodic traveling waves. 

e. Spontaneous sources as a precursor of spontaneous 
periodic waves The periodic standing waves emerging 
via the cross-diffusion driven instability described above, 
are in turn unstable themselves. Fig. 4(d) shows a con- 
tinuation of simulation of fig. 4(c). The standing waves 
are observed for a long time, as they are stable within 
the space of functions with spatial period 2ir/k = L/6, 
and the numerical initial conditions are almost exactly 
periodic with that period, up to small errors resulting 
from finite precision arithmetics. The small symmetry- 
breaking numerical errors allow for an instability of the 
periodic standing waves to develop, during which some of 
the standing waves occur later than others. When this in- 
stability sufficiently develops, there is a sudden, "hard" 
transition to propagating waves. The spatial period of 
the propagating waves is twice longer than the spatial 
period of preceding standing waves. We stress that the 
traveling waves do not appear via anything like "bifurca- 
tion" from standing waves, at least in the examples we 
considered. 

Notice that the long transient solution shown in 
fig. 4(c,d) is a periodic standing wave by its symmetry, 
but it also looks like a periodic set of focal sources, syn- 
chronously sending out solitary waves which then anni- 
hilate each other. As can be seen in fig. 1, apart from 
the symmetry, this sort of transient before the onset of 
periodic waves is typical, and only its duration varies in 
different simulations. That is, the special initial condi- 
tions in fig. 4(c,d) only affect the symmetry and the du- 
ration of the transient, but not its qualitative character. 
A similar route to traveling waves via unstable periodic 
set of "focal sources" standing waves is obseved in the 
RM model. 



V. CONCLUSION 

The considered examples demonstrate an unusual type 
of behaviour. The systems are oscillatory, but the spa- 
tially uniform oscillations are unstable. The systems can 
also demonstrate standing periodic waves, which are also 



unstable. These instabilities lead to periodic propagating 
waves, which seems to be the only stable regime. This 
regime emerges spontaneously even when boundary con- 
ditions disallow propagating waves. The periods of the 
waves can be in a certain interval with strict boundaries, 
both upper and lower. Nearer the upper end of the in- 
terval, i.e. at longer wavelengths, the periodic waves do 
not propagate steadily but are modulated. Transition 
from steady to modulated propagation is soft and has 
empirical features of a supercritical Hopf bifurcation (of 
a relative equilibrium), i.e. possibly an Eckhaus mecha- 
nism. 

The defining features described above are sufficiently 
generic, and the phenomenon of spontaneous periodic 
traveling waves does not disappear as the parameters are 
varied, nor it is restricted just to one model. This be- 
haviour does not fall into existing classification of pattern 
formation scenarios. The detailed mechanisms of spon- 
taneous generation and maintenance of periodic travel- 
ing waves require further investigation. However, it is 
clear that cross-diffusion is an essential factor, since its 
replacement with, or adding of significant amount of, self 
diffusion eliminates the effect. Cross-diffusion phenom- 
ena are known in a in a variety of physical situations. 
For example, spontaneous periodic waves have been ob- 
served in a Burridge-Knopoff mathematical model of 
earthquakes [24, 25]. That model belongs to the class (1), 
with only one nonzero element of matrix D, as in our sim- 
ulations shown in fig. 1(a) and (g) but constant, and ex- 
citable FitzHugh-Nagumo local kinetics. It is not known 
whether the spontaneous waves in the Burridge-Knopoff 
model have a finite interval of allowed wavenumbers, as 
illustrated by fig. 3 for our case, however other described 
features of those waves are similar to those described here 
and are likely to have a similar nature. Further investi- 
gation of the mechanism of generation of such waves is 
a subject for further study which is of broad physical in- 
terest as a new pattern forming mechanism in dissipative 
spatially distributed systems. 

Returning to the application that originally motivated 
this study, attempts to explain waves observed in cyclic 
biological populations, using reaction-diffusion models, 
had to involve spatially- nonuniform external factors, e.g. 
sites of increased mortality due to environmental condi- 
tions [14, 26]. Such factors are needed to disallow uni- 
form oscillations. Our present results imply that such 
factors may not be necessary if cross-diffusion interac- 
tion are taken into account as the uniform oscillations 
may be unstable and waves form spontaneously. 
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